!     Copyright (C) 2010 Benjamin Piaud
!
!     LIMBES is free software; you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation; either version 3, or (at your option)
!     any later version.

!     LIMBES is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.

!     You should have received a copy of the GNU General Public License
!     along with LIMBES; if not, see <http://www.gnu.org/licenses/>


module LIMBES_mesh_mod
  integer*4 :: nx,ny,Nmax,fluid_ntot,b1_ntot, b2_ntot
  integer*4 :: VTK_xmin,VTK_xmax,VTK_ymin,VTK_ymax
  integer*4 ,allocatable :: phase(:,:),wall_RLUD(:,:)
  integer*4 ,allocatable :: cell_xi(:),cell_yi(:)
  integer*4 ,allocatable:: cell_ni(:,:),neighbour_RLUD(:,:)
  real*8 ,allocatable :: alpha1(:,:),alpha2(:,:) 
  character*2, allocatable :: boundary_type1(:), boundary_type2(:,:)
  real*8 :: delta_x,zero
  parameter(zero=0.0D+0)

contains

  subroutine boundary()
    implicit none
    integer*4 :: xi,yi,ni,diffL,diffR,diffD,diffU,nb1,nb2,wall

    allocate(wall_RLUD(1:fluid_ntot,1:4))
    !     detecting boundaries nodes
    !     wall_RLUD(ni,1)=1 means a wall at the Right side of cell ni
    !     wall_RLUD(ni,2)=1 means a wall at the Left side of cell ni
    !     wall_RLUD(ni,3)=1 means a wall at the Up side of cell ni
    !     wall_RLUD(ni,4)=1 means a wall at the Down side of cell ni
    do xi=1,nx
       do yi=1,ny
          if (phase(xi,yi).EQ.1) then
             ni=cell_ni(xi,yi)
             if (xi.EQ.1) then
                diffL=0
             else
                diffL=iabs(phase(xi,yi)-phase(xi-1,yi))
             endif
             if (xi.EQ.nx) then
                diffR=0
             else
                diffR=iabs(phase(xi+1,yi)-phase(xi,yi))
             endif
             if (yi.EQ.1) then
                diffD=0
             else
                diffD=iabs(phase(xi,yi)-phase(xi,yi-1))
             endif
             if (yi.EQ.ny) then
                diffU=0
             else
                diffU=iabs(phase(xi,yi+1)-phase(xi,yi))
             endif
             wall_RLUD(ni,1)=0
             wall_RLUD(ni,2)=0
             wall_RLUD(ni,3)=0
             wall_RLUD(ni,4)=0
             if (diffL.EQ.1) wall_RLUD(ni,2)=1   
             if (diffR.EQ.1) wall_RLUD(ni,1)=1
             if (diffD.EQ.1) wall_RLUD(ni,4)=1
             if (diffU.EQ.1) wall_RLUD(ni,3)=1
          endif
       enddo
    enddo

    !     count boundary nodes with 1 wall and 2 walls
    nb1=0
    nb2=0
    do ni=1,fluid_ntot
       wall=wall_RLUD(ni,1)+ wall_RLUD(ni,2) + wall_RLUD(ni,3)+ wall_RLUD(ni,4)
       if (wall.GT.2) then
          write(*,*) 'bad boudary entry (only 2 walls/cell allowed)'
          stop
       endif
       if (wall.EQ.1) nb1=nb1+1
       if (wall.EQ.2) nb2=nb2+1
    enddo
    b1_ntot=nb1
    b2_ntot=nb2
    
    allocate(boundary_type1(1:fluid_ntot))
    allocate(alpha1(1:fluid_ntot,1:4))
    allocate(boundary_type2(1:fluid_ntot,1:2))
    allocate(alpha2(1:fluid_ntot,1:4))
    
  end subroutine boundary


  subroutine neighbour()
    implicit none
    integer*4 :: ni,xi,yi,x1,y1,n1
    
    allocate(neighbour_RLUD(1:fluid_ntot,1:4))
    !     evaluating node neighbours in each axes 
    !     (including periodic boundary condition)
    !     neighbour_RLUD(ni,1) Right neigbour of node ni
    !     neighbour_RLUD(ni,2) Left neigbour of node ni   
    !     neighbour_RLUD(ni,3) Up neigbour of node ni
    !     neighbour_RLUD(ni,4) Down neigbour of node ni      
    do xi=1,nx
       do yi=1,ny
          if (phase(xi,yi).EQ.1) then
             ni=cell_ni(xi,yi)
             !     Right Neighbour
             x1=xi+1
             if (x1.GT.nx) x1=mod(x1,nx) !periodic boundary condition
             n1=cell_ni(x1,yi)
             if (wall_RLUD(ni,1).EQ.1) n1=ni
             neighbour_RLUD(ni,1)=n1
             !     Left Neighbour
             x1=xi-1
             if (x1.LT.1) x1=mod(x1,nx)+nx !periodic boundary condition
             n1=cell_ni(x1,yi)
             if (wall_RLUD(ni,2).EQ.1) n1=ni
             neighbour_RLUD(ni,2)=n1               
             !     Up Neighbour
             y1=yi+1
             if (y1.GT.ny) y1=mod(y1,ny) !periodic boundary condition
             n1=cell_ni(xi,y1)
             if (wall_RLUD(ni,3).EQ.1) n1=ni
             neighbour_RLUD(ni,3)=n1
             !     Down Neighbour
             y1=yi-1
             if (y1.LT.1) y1=mod(y1,ny)+ny !periodic boundary condition
             n1=cell_ni(xi,y1)
             if (wall_RLUD(ni,4).EQ.1) n1=ni
             neighbour_RLUD(ni,4)=n1
          endif
       enddo
    enddo

  end subroutine neighbour

  subroutine write_mesh
    implicit none
    integer*4 :: ni
    !     writing the mesh information
    open(10,file='./mesh.in')
    write(10,*) '#/number of cells/ ','/delta x/','/nx/','/ny/'
    write(10,*) '# /xmin/  /xmax/  /ymin/   /ymax/  for VTK output'
    write(10,*) '#/cell index/ ','/xi cell/ ','/yi cell/ ','/cell on Right/ ','/cell on Left/ ','/cell on Up/ ','/cell on Down/'   
    write(10,*) fluid_ntot, delta_x, nx, ny
    write(10,*) VTK_xmin,VTK_xmax,VTK_ymin,VTK_ymax
    do ni=1,fluid_ntot
       write(10,'(11i7)') ni,cell_xi(ni),cell_yi(ni),neighbour_RLUD(ni,1),neighbour_RLUD(ni,2),neighbour_RLUD(ni,3),neighbour_RLUD(ni,4)
    enddo
    close(10)
  end subroutine write_mesh
  
  subroutine write_boundary
    implicit none
    integer*4 :: nb1,nb2,ni,wall
    !     writing boundary information
  nb1=0
  nb2=0
  open(11,file='./boundary.in')
  write(11,*) '#number of cells with 1 wall /','#number of cells with 2 walls'
  write(11,*) '#/boundary index/ ','/boundary type/ ','/cell index/ ','/Right wall/ ', '/Left wall/ ','/Up wall/ ', '/Down wall/ ','/alpha/' 
  write(11,*) b1_ntot,b2_ntot
  do ni=1,fluid_ntot
     wall=wall_RLUD(ni,1)+ wall_RLUD(ni,2)+ wall_RLUD(ni,3)+ wall_RLUD(ni,4)
     if (wall.EQ.1) then
        nb1=nb1+1
        write(11,*) nb1,'"',boundary_type1(ni),'"',ni,wall_RLUD(ni,1), wall_RLUD(ni,2),wall_RLUD(ni,3),wall_RLUD(ni,4),&
             alpha1(ni,1),alpha1(ni,2),alpha1(ni,3),alpha1(ni,4)
     endif
  enddo
  write(11,*) '#/boundary index/ ','/boundary type 1/ ','/boundary type 2/ ','/cell index/ ','/Right wall/ ', '/Left wall/ ','/Up wall/ ', '/Down wall/ ','/alpha/'
  do ni=1,fluid_ntot
     wall=wall_RLUD(ni,1)+ wall_RLUD(ni,2) + wall_RLUD(ni,3)+ wall_RLUD(ni,4)
     if (wall.EQ.2) then
        nb2=nb2+1
        write(11,*) nb2,'"',boundary_type2(ni,1),'" ','"',boundary_type2(ni,2),'"',ni, &
             wall_RLUD(ni,1),wall_RLUD(ni,2),wall_RLUD(ni,3),wall_RLUD(ni,4),&
             alpha1(ni,1),alpha1(ni,2),alpha1(ni,3),alpha1(ni,4),&
             alpha2(ni,1),alpha2(ni,2),alpha2(ni,3),alpha2(ni,4)
     endif
  enddo
  close(11)
end subroutine write_boundary

end module LIMBES_mesh_mod
